Large-scale simulations of fluctuating biological membranes 
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We present a simple, and physically motivated, coarse-grained model of a lipid bilayer, suited 
for micron scale computer simulations. Each ^ 25nm^ patch of bilayer is represented by a spher- 
ical particle. Mimicking forces of hydrophobic association, multi-particle interactions suppress the 
exposure of each sphere's equator to its implicit solvent surroundings. The requirement of high 
equatorial density stabilizes two-dimensional structures without necessitating crystalline order, al- 
lowing us to match both the elasticity and fluidity of natural lipid membranes. We illustrate the 
model's versatility and realism by characterizing a membrane's response to a prodding nanorod. 



INTRODUCTION 

Lipid bilayers form the basis of biological membranes. 
Integral membrane proteins and a variety of small 
molecules are embedded in this two-dimensional fluid, 
which is stabilized by hydrophobic interactions: the bi- 
layer structure effectively shields the lipid's hydropho- 
bic alkane chains from exposure to the aqueous solvent. 
Despite this complexity on the molecular level, biologi- 
cal membranes at large length scales are well character- 
ized by surprisingly few material properties. In partic- 
ular, their behavior on large length scales is consistent 
with that of crude elastic models, which take as input 
only a bending rigidity, i.e., the membrane's resistance 
to smooth shape deformations. It is unclear below what 
scale such representations become inappropriate, in part 
due to computational difficulties of incorporating ther- 
mal fluctuations, which invariably gain importance as the 
scale of observation is reduced from the macroscopic. At 
the opposite extreme computer simulations of lipid bi- 
layers at atomistic resolution provide detailed access to 
the spectrum of thermal fluctuations. They are limited, 
however, to length and time scales of tens of nanometers 
and hundreds of nanoseconds, respectively 

Many important biological phenomena involve mem- 
brane deformations over several micrometers and span 
seconds or even minutes. They thus occur on length 
scales intermediate between those natural to elastic con- 
tinuum models and to atomistic representations. In at- 
tempts to bridge this gap, a large variety of simplified 
models for interactions between lipid molecules have been 
proposed in the literature. These range from systemati- 
cally coarse-grained systems [l, Q to solvent-free heuris- 
tic models 0, HI ; a recent review can be found in Ref. @. 
Common to these models is an attempt to mimic the am- 
phipathic character of individual lipid molecules. While 
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considerably less expensive than atomistic simulations in 
their computational demands, these approaches are still 
limited in the scope of fluctuations and response they can 
feasibly capture. 

Extending computer simulations to examine large scale 
behaviors such as aggregation of membrane-associated 
proteins or deformations induced by growth of an actin 
network would appear to require coarse-graining beyond 
the scale of individual lipid molecules. Several such ap- 
proaches have been proposed 0, 0, 0]- Most represent 
in a discrete way the fluctuations implied by Helfrich's 
continuum model of an elastic sheet [13, [HI • These re- 
quire estimating local curvature of a discretized mani- 
fold, which can be both numerically unstable and taxing 
to implement. Furthermore, these approaches do not at- 
tempt to make a close connection with the physical forces 
responsible for membrane cohesion and elasticity. 

In this work we propose a new simulation model at the 
many-lipid scale that follows transparently from the sta- 
tistical mechanics of hydrophobicity and the underlying 
membrane thermodynamics. Despite its simplicity, our 
model successfully captures several important properties 
of lipid bilayer physics, such as intrinsic fluidity and the 
ability to spontaneously assemble into two-dimensional 
sheets. It is furthermore sufficiently versatile to repro- 
duce a wide range of biologically relevant elastic proper- 
ties. 

In our model, we envision the lipid bilayer as a collec- 
tion of small membrane patches of size d ^ 5 nm, roughly 
the thickness of a typical bilayer [lH, each comprising 
~ 100 phopholipid molecules. For geometric simplicity, 
we represent each patch as a volume-excluding sphere 
with an axis of rotational symmetry pointing from one 
polar head group region to the other. Cohesion of such 
patches is due of course to the presence of water: Expos- 
ing the hydrophobic portion (i.e., the equatorial region of 
our model spheres) to solvent incurs a free energetic cost, 
while exposing the hydrophilic portion (i.e., the polar 
caps of our model spheres) is thermodynamically advan- 
tageous. At length scales ^ 1 nm both of these contri- 
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butions should be proportional to the exposed area [l3| . 
Remarkably, these considerations alone are sufficient to 
successfully mimic the flexibility and fluidity of natural 
bilayers. 

The model we have developed from these simple phys- 
ical notions resembles in some respects one reported long 
ago by Leibler and coworkers [ij]. They similarly con- 
sidered association of spherical units each representing a 
bilayer patch comparable in size to the membrane's thick- 
ness. The anisotropic interactions acting among their 
particles, however, were devised not to reflect pertinent 
thermodynamic driving forces at this length scale, but 
instead to foster formation of fluid elastic sheets. In 
our view the potential energy function used in that work 
would be difficult to motivate on microscopic grounds, 
and no attempt was made in Ref. Q to do so. For this 
reason we expect that our approach will generalize more 
naturally to describe scenarios that involve membrane 
properties beyond long-wavelength fluidity and elastic- 
ity. 

As a more immediate and practical justiflcation for our 
new approach, we found that the original implementation 
of the model, as described in Ref. [3, does not yield stable 
two-dimensional structures. In the Appendix we detail 
a revised version of that model which is consistent with 
previously reported properties. But even in this case the 
model membrane's bending rigidity is atypically small for 
biophysical systems. 



MODEL 

For a collection of N particles, each representing a 
patch of lipid bilayer, we adopt the energy function 



N 



(1) 



where the hard-core potential Uuc enforces the con- 
straint that the separation between any two particles is at 
least d. The quantities neq and n^'^j characterize, respec- 
tively, the equatorial and polar coordination numbers of 
particle i. The functions Aeq(n) and Apo\{n) determine 
solvent exposure of these two regions based on their co- 
ordination. The positive constant e sets the scale of these 
solvent-mediated interactions. Based on the surface ten- 
sion between water and oil, 7 ~ 50mJ/m^, and the hy- 
drophobic surface area of a membrane patch, A ^ 60nm^, 
we expect ^ 740/cbT to be an appropriate value for 
e. 

In detail, we define the fluctuating density of a parti- 
cle's equatorial and polar neighbors as 
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The contributions of particle j to the coordination densi- 
ties of particle i are thus determined both by the distance 
rij = \Tij\ between their centers, and by the normalized 
projection Zij = Vij • d^/r^j of their separation vector Vij 
onto the axis of particle i (which points along the unit 
vector di). They are attenuated by the functions 



Geq(r) Gpoi(r) < 



if r < Ta, 
if Ta < r < 
otherwise, 
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and 



ffeq(^') = l-ffpol(^') 
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otherwise. 
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over scales determined by parameters r^j, r^, z^^ and Zb. 
In this way, particle j counts toward the equatorial (po- 
lar) coverage of particle i only if it lies near its center and 
not too far off the equatorial plane (polar axis). 

For a partially occluded object, such as our idealized 
membrane particles when surrounded by neighbors, cal- 
culating the area accessible to solvent molecules of finite 
size is a nontrivial operation. Furthermore, this task is 
not sensible to carry out in detail for the reduced rep- 
resentation we have chosen. Essential features of the 
functions Aeq(n) and Apoi{n) are nonetheless straightfor- 
ward to ascertain: As n increases, exposed area at first 
declines steadily. At some value n lower than the max- 
imum n("^^^) permitted by steric constraints, it should 
nearly vanish, since perfect close-packing is not needed 
to thoroughly exclude solvent from the bilayer's interior 
(or from the polar exterior). For n > n, variation in 
exposed area should be very weak. We caricature this 
dependence with computationally inexpensive, piecewise 
linear functions: 



Aeq{n) = 1 - — (n < neq), (n > neq), 

^eq 



(6) 



^poi(^) = 1 - - — {n < npoi), (n > npoi). (7) 

Appropriate values of neq and npoi will depend on the rel- 
ative thicknesses of hydrophobic and polar regions, spec- 
ified in our model by Za and Z5. 

The membrane free energy we have described is 
markedly multi-body in character, but in a way that is 
both physically meaningful and simple to understand. Its 
parameters correspond intuitively to the geometry and 
chemistry of constituent molecules. Below we present re- 
sults of Monte Carlo computer simulations for this model, 
which employed single particle translations and rotations 
as trial moves, as well as shape deformations of the peri- 
odically replicated simulation box when an external ten- 
sion was imposed. For these calculations we selected 
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zl 




0.01 


0.2 


126.8±0.8 


0.02 


0.2 


63.4±0.5 


0.05 


0.2 


25.2±0.1 


0.1 


0.2 


12.1±0.1 


0.4 


0.6 


2.1±0.03 



TABLE I: Computed bending rigidities for different values of 
the model parameters and zl. 
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FIG. 1: (a) Illustration of our model. Each particle corre- 
sponds to a fragment of lipid bilayer, comprised of a cen- 
tral hydrophobic core and two surrounding hydrophilic layers. 
The unit vector specifies the orientation of particle i. Also 
shown are the spatial regions used to compute the numbers 
of equatorial (enclosed by dashed line) and polar (enclosed by 
dotted line) neighbors of particle i. (b) Random initial con- 
figuration for a trajectory of = 864 particles. As time pro- 
gresses, particles quickly form two-dimensional patches that 
continue to coarsen. Shown are snapshots after 100 (c) and 
1000 (d) Monte Carlo sweeps. 
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Va/d = l.?>,n/d = 1.7, zl = 0.05, zl = 0.2, n, 
np = 1, which yield a rigidity typical of biological mem- 
branes. By varying these values, it is possible to tune the 
elasticity of the assembled bilayer, as well as more subtle 
properties like internal viscosity or rupture tension. 



RESULTS 

We first demonstrate that a sheet-like configuration 
is indeed the equilibrium state of our model at finite 
temperature. Fig. [T] shows snapshots from a Monte 
Carlo trajectory in which an initially disperse collection 
of membrane particles spontaneously organizes into two- 
dimensional structures that then diffuse and coalesce [22[ . 
This coarsening process is expected to proceed until only 
a single membrane sheet remains in the simulation box. 
A free boundary in the resulting structure can be avoided 
either by forming a membrane sheet that spans the pe- 
riodically replicated simulation box, or by adopting a 
boundary-free geometry such as a spherical vesicle. 

We next provide evidence that the elastic properties 
of our model match quantitatively those of natural lipid 
bilayer s. Specifically, on length scales well beyond a par- 



ticle radius d (corresponding to the bilayer 's thickness), 
its shape fluctuations are well described by the Helfrich 
model of incompressible fluid elastic sheets [lol, [HI , with 
macroscopic material properties in accord with experi- 
mental measurements. In Helfrich's model, a nearly flat 
segment of membrane exhibits a fluctuation spectrum 
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Here, hq = h{x) exp {—iq ■ x) dx is the Fourier trans- 
form of the membrane height h{x), gr is a two-dimensional 
wavevector conjugate to the Cartesian position x on a 
flat reference membrane with area A, and (.) denotes 
the equilibrium average. For lipid mixtures common in 
cell membranes, the bending rigidity ranges from 10 
to 30 /cb^ [15]. To the extent that the membrane's area 
per molecule is constant, the surface tension a plays the 
role of minus chemical potential and increases roughly 
linearly with applied lateral tension r. 

For the purpose of estimating the normal mode fluctua- 
tions of Eq.[8]from simulation, we wish to avoid imposing 
lateral tension (requiring r = 0). As a practical matter, 
however, it is convenient in this calculation to prescribe 
a fixed box geometry (and thus a fixed set of wavevec- 
tors). We began by performing simulations in which the 
box size was allowed to fluctuate at zero lateral tension. 
The resulting average box size was then adopted as a 
constraint for a second set of simulations, in which we 
computed the distribution of height fluctuations [23^ . At 
large length scales we indeed observe the characteristic 



he 



(X q ^ behavior predicted by the Helfrich model. 



as shown in Fig. O The computed proportionality coef- 
ficient indicates a bending rigidity of ~ 25/cbT, well 
within the range of measured values ^[ . By manipulat- 
ing the parameters z^ and zl, we are able to reproduce 
the elastic behavior of membrane sheets over a wide range 
of bending rigidities (Table [Tj). 

The advantage of our model lies in a facile ability to 
address mesoscale response without sacrificing the mi- 
croscopic basis of corresponding fluctuations. As a rep- 
resentative biophysical example that calls for these ca- 
pabilities, we considered the resistance of a fluctuating 
membrane to impingement of a nanorod oriented perpen- 
dicular to the lipid bilayer. Experimental realizations of 
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FIG. 2: Spectrum of height fluctuations around a flat refer- 
ence state, equilibrated at zero lateral tension. The dashed 
line is a fit to the expected behavior (|8| with a = 0, yielding 
a bending rigidity n — 25.2/cbT. 



this situation include extension of polymerizing actin fil- 
aments close to a cell membrane [16| and external forcing 
of a carbon nanotube against a cell wall 



Insets of Fig. [3] depict the reversible membrane defor- 
mations we have studied in this context. The prodding 
nanorod is modeled here as a volume-excluding, rigid 
spherocylinder of radius R = 3d, which defines a region 
in space inaccessible to the membrane particles. Its ver- 
tical displacement / determines the size of the membrane 
deformation. We set / = for a nanorod that would con- 
tact the membrane, in a completely flat configuration, at 
a single point. To prevent global translation of the mem- 
brane when / > 0, we constrain the vertical positions of a 
small number of membrane particles. This pinning could 
be viewed as a mimicry of cytoskeletal attachments that 
would suppress overall translations in a living cell [l2[ . 

To calculate membrane-induced forces on the nanorod, 
we treat the rod height / as a dynamical variable subject 
to an external potential, V{1) = {1/2)K {I — /q)^, in addi- 
tion to the fluctuating constraints imposed by membrane 
particles. For several values of Iq we computed the force 
on a protrusion of length I from simulations as 



{I) - lo 



{5iy ) - ksT/K 



{siy 



((0-0 



, (9) 



which follows from an expansion of the membrane free 
energy to quadratic order in / around I. 

Force-extension curves for this extrusion process are 
plotted in Fig. [3] for three different values of r. The 
restoring force initially increases in proportion to filament 
length, and ultimately reaches a plateau value at large I. 
This limiting force is in good agreement with the con- 
stant force 2i\\J2kg predicted by the Helfrich model for 
the stretching deformation of an elastic cylinder [11, [l^, 
provided that a is dominated by lateral tension. Simi- 
lar agreement is found for the diameter of the membrane 




FIG. 3: Restoring force as a function of protrusion length 
for three different values of the lateral tension r: 0.0002e/(i^ 
(squares), 0.0005e/(i^ (circles), and O.OOle/d^ (triangles). For 
each value of r and /o we performed computer simulations at 
two values of K (0.2 and 1 in units of e/d^). The force was 
evaluated using (|9} at a length / corresponding to the arith- 
metic mean of the two average lengths (/) . The solid lines 
show the result 2i^\/2kg expected for long cylindrical protru- 
sions in the limit of zero temperature with g — r [isl . Fl9j |. 
The inset shows cutaway views of typical configurations at 
T = O.OOOSe/d^ and h = 17.5,29.5,41.5 (left to right), illus- 
trating the transition from a global deformation to a localized, 
tubular protrusion. In our calculations the nanorod is placed 
at the center of the periodically replicated simulation box, 
and a small number of membrane particles are immobilized 
at the boundary to avoid overall membrane translation. 



tubule (data not shown). At intermediate extensions, 
/ exhibits a local maximum, signaling a mechanical in- 
stability associated with the transformation between two 
different classes of membrane configurations. This ob- 
servation is consistent with recent experiments on the 
formation of membrane tethers from giant vesicles l2Cl 
and conclusions drawn from continuum calculations [l8| . 



OUTLOOK 

In addition to mechanical responses exemplified here 
by nanorod protrusion, several other biologically relevant 
membrane processes could addressed through straight- 
forward extensions of our model. For example, lateral 
demixing of multicomponent lipid bilayers could be stud- 
ied by endowing each membrane particle with a degree of 
freedom representing the local composition. The organi- 
zation of embedded proteins could be simulated using a 
description of the macromolecule that is compatible with 
our membrane model. Studying time-dependent behav- 
ior would require in addition a set of dynamical rules 
for updating particle arrangements. Most simply, the 
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Monte Carlo trajectories we have used here to sample 
equilibrium configurations could be interpreted as an ap- 
proximation of physical time evolution. Integrating equa- 
tions of motion that involve derivatives of potential en- 
ergy would require smoothed versions of the interactions 
we have described. Already, the mesoscale realism and 
computational economy of our approach recommends its 
use for examining a wide range of lipid bilayer fiuctuation 
phenomena beyond the reach of molecular models. 

We thank Anthony Maggs for helpful discussions of the 
work presented in Ref. [ij. This work was supported by 
the Director, Office of Science, Office of Basic Energy Sci- 
ences, of the U.S. Department of Energy under Contract 
No. DE-AC02-05CH11231. 



APPENDIX 

The model proposed by Drouffe and coworkers in 
Ref. [3 operates on the same lengthscale and focuses on 
the same degrees of freedom as the model presented in 
this work. It employs an energy functional of the form 



U — UuC + ^an + ^den- 



(10) 



The first term is a pairwise repulsive energy that penal- 
izes overlap of different particles, which we treat as hard 
spheres. The second term is a pairwise anisotropic inter- 
action energy. 
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(11) 

where Bivij) is a positive weight function that limits the 
range of the interaction to ^ 2d^ r] = 1 in Ref. Q, ol G 
{1,2} depending whether the particles are symmetric, 
and 



g{x) = 0.75x^ + 0.25X + 0.8 



(12) 



is a monotonically increasing function over the range of 
its argument x G [0, 1]. 

The last term in (p!Q|) is a multibody potential that 
favors configurations in which each particle is surrounded 
by 6 nearest neighbors. 



^den = e^C{pi), 



(13) 



where C{pi) = {pi - 6)^, pi = Y^j^iKvij), and h{rij) 
is a weight function with unit value for small particle 
separations. 

To implement this model one must specify the depen- 
dence of B{r) and h{r) on inter-particle distance, which 
was done only graphically in Ref. We assign these 



quantities simple functional forms that are piecewise lin- 
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ear m r 



B{r) 



h{r) 



1.5, 



if r < 1.85, 



^2^gL°i.85.\ if 1.85 < r < 2.0, 

0, otherwise, 

1, ifr<1.6, 
^05^^ if 1.6 < r < 1.95, 
0, otherwise. 



(14) 



(15) 



Values of B and h resulting from this choice closely re- 
semble those plotted in Ref. 14. 

We performed Monte Carlo computer simulations of 
the model defined by (p!Q |) - (p!5|) for a = 1 at a tem- 
perature T = 1.6e. For these conditions, Drouffe and 
coworkers report the spontaneous assembly of randomly 
distributed and oriented particles into a two-dimensional 
sheet, in which the orientation vectors of neighboring 
particles are nearly parallel. However, in our simula- 
tions we find that for these parameter values the system 
forms many small clusters of particles without apparent 
internal order. Similar configurations are obtained when 
a completely ordered sheet is used as the initial condi- 
tion (Fig. [4^), which indicates that the inability to form 
stable membranes is thermodynamic in origin, and not 
due to dynamical artefacts caused by the Monte Carlo 
method. 

The lack of stable membrane-like structure in this pub- 
lished version of the model is straightforward both to un- 
derstand and to remedy. When 77 = 1, the first term in 
the brackets of Eq. ([TT]) actually disfavors parallel align- 
ment of neighboring particles' orientational vectors; this 
energy contribution is minimum when axes of adjacent 
particles are antiparallel (in the case a = 1) or perpen- 
dicular {a = 2). Suspecting a typographical error in 
Ref. [3, we set 77 = — 1 in order to instead favor the de- 
sired alignment. Equilibrium states resulting from this 
modification, however, still feature a collection of many 
small aggregates of p articles at the temperature T = 1.6 
reported in Ref. |lj. In this case the adhesive interac- 
tion (fT3|) is insufficient to overcome the loss in entropy 
associated with forming a single membrane sheet. Only 
after reducing the temperature to T = 0.5e were we able 
to observe stable sheet-like structures and vesicles in our 
simulations (Fig. [IJd). This modified set of model param- 
eters results in membrane sheets with a bending rigidity 
tv = 4.9±0.1/cbT, which is smaller than values measured 
for biophysical phospholipid bilayers ^Ji] . 
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